{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bisect import bisect_right\n",
    "\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 插值\n",
    "## 拉格朗日插值"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10.722755505364201\n"
     ]
    }
   ],
   "source": [
    "def lagrange_interploate(xs, ys):\n",
    "    \"\"\"\n",
    "    拉格朗日插值\n",
    "    :param xs: 插值节点 x\n",
    "    :param ys: 插值节点 y\n",
    "    :return: 插值多项式\n",
    "    \"\"\"\n",
    "    def do_interploate(x):\n",
    "        ret = 0\n",
    "        for i in range(len(ys)):\n",
    "            l = 1\n",
    "            for j in range(len(xs)):\n",
    "                if j == i:\n",
    "                    continue  # 注意这句\n",
    "                l *= (x - xs[j]) / (xs[i] - xs[j])\n",
    "            ret += l * ys[i]\n",
    "\n",
    "        return ret\n",
    "\n",
    "    return do_interploate\n",
    "\n",
    "\n",
    "xs = [100, 121, 144]\n",
    "ys = [10, 11, 12]\n",
    "func = lagrange_interploate(xs, ys)\n",
    "print(func(115))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 牛顿插值"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8.555555555555557\n"
     ]
    }
   ],
   "source": [
    "def newton_interploate(xs, ys):\n",
    "    \"\"\"\n",
    "    牛顿插值\n",
    "    :param xs: 插值节点 x\n",
    "    :param ys: 插值节点 y\n",
    "    :return: 插值多项式\n",
    "    \"\"\"\n",
    "    dq_tab = np.zeros((len(xs), len(xs)))  # DP 计算差商表\n",
    "    dq_tab[0] = ys\n",
    "    for i in range(1, len(xs)):\n",
    "        for j in range(len(xs) - i):\n",
    "            dq_tab[i][j] = (dq_tab[i - 1][j] - dq_tab[i - 1]\n",
    "                            [1 + j]) / (xs[j] - xs[j + i])\n",
    "\n",
    "    def do_interploate(x):\n",
    "        x_mul = 1  # 计算 (x-x0)(x-x1)...(x-x_{i-1})\n",
    "        ret = 0\n",
    "        for i in range(len(xs)):\n",
    "            ret += dq_tab[i][0] * x_mul\n",
    "            x_mul *= (x - xs[i])\n",
    "\n",
    "        return ret\n",
    "\n",
    "    return do_interploate\n",
    "\n",
    "\n",
    "xs = [0, 1, 4, 3, 6]\n",
    "ys = [0, -7, 8, 5, 14]\n",
    "func = newton_interploate(xs, ys)\n",
    "print(func(5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hermite 插值"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-4.0\n"
     ]
    }
   ],
   "source": [
    "def hermite_interploate(xs, ys, ys2):\n",
    "    \"\"\"\n",
    "    Hermite 插值\n",
    "    :param xs: 插值节点 x\n",
    "    :param ys: 插值节点 y\n",
    "    :param ys2: 插值节点 y 导数\n",
    "    :return: 插值多项式\n",
    "    \"\"\"\n",
    "    def do_interpolate(x):\n",
    "        ret = 0\n",
    "        for i in range(len(xs)):\n",
    "            l = 1  # l_i(x)\n",
    "            l2 = 0  # l_i(x_i) 的导数\n",
    "            for j in range(len(xs)):\n",
    "                if j == i:\n",
    "                    continue\n",
    "                l2 += 1 / (xs[i] - xs[j])\n",
    "                l *= (x - xs[j]) / (xs[i] - xs[j])\n",
    "            alpha = (1 - 2 * (x - xs[i]) * l2) * l ** 2\n",
    "            beta = l ** 2 * (x - xs[i])\n",
    "            ret += ys2[i] * beta + ys[i] * alpha\n",
    "\n",
    "        return ret\n",
    "\n",
    "    return do_interpolate\n",
    "\n",
    "\n",
    "xs = [1, 2]\n",
    "ys = [2, 3]\n",
    "ys2 = [1, -1]\n",
    "func = hermite_interploate(xs, ys, ys2)\n",
    "print(func(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 分段线性插值"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def frac_linear_interpolate(xs, ys):\n",
    "    \"\"\"\n",
    "    分段线性插值\n",
    "    :param xs: 插值节点 x\n",
    "    :param ys: 插值节点 y\n",
    "    :return: 插值多项式\n",
    "    \"\"\"\n",
    "    def do_interpolate(x):\n",
    "        ind = bisect_right(xs, x)  # bisect 的使用\n",
    "        return ys[ind] * (x - xs[1 + ind]) / (xs[ind] - xs[1 + ind]) + ys[1 + ind] * (x - xs[ind]) / (xs[1 + ind] - xs[ind])\n",
    "\n",
    "    return do_interpolate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 线性拟合"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.381081081081081\n"
     ]
    }
   ],
   "source": [
    "def linear_regression(xs, ys):\n",
    "    \"\"\"\n",
    "    线性拟合\n",
    "    :param xs: 拟合数据 x\n",
    "    :param ys: 拟合数据 y\n",
    "    :return: 拟合结果\n",
    "    \"\"\"\n",
    "    n = len(xs)\n",
    "    xy_sum = np.sum(xs * ys)  # np 的使用\n",
    "    x_sum = np.sum(xs)\n",
    "    y_sum = np.sum(ys)\n",
    "    x2_sum = np.sum(xs ** 2)\n",
    "\n",
    "    b = (n * xy_sum - x_sum * y_sum) / (n * x2_sum - x_sum ** 2)\n",
    "    a = (y_sum - b * x_sum) / n\n",
    "\n",
    "    def do_regression(x):\n",
    "        return x * b + a\n",
    "\n",
    "    return do_regression\n",
    "\n",
    "\n",
    "xs = np.array([0, 1, 2, 3, 5])\n",
    "ys = np.array([1.1, 1.9, 3.1, 3.9, 4.9])\n",
    "func = linear_regression(xs, ys)\n",
    "print(func(4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 数值积分"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 机械求积公式"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.75\n",
      "0.8353317053317053\n",
      "0.8425925925925926\n"
     ]
    }
   ],
   "source": [
    "def T(f, a, b):\n",
    "    \"\"\"\n",
    "    梯形公式\n",
    "    :param f: 被积函数\n",
    "    :param a: 积分下限\n",
    "    :param b: 积分上限\n",
    "    :return: 积分结果\n",
    "    \"\"\"\n",
    "    return (b - a) * (f(a) + f(b)) / 2\n",
    "\n",
    "\n",
    "def S(f, a, b):\n",
    "    \"\"\"\n",
    "    辛普森公式\n",
    "    :param f: 被积函数\n",
    "    :param a: 积分下限\n",
    "    :param b: 积分上限\n",
    "    :return: 积分结果\n",
    "    \"\"\"\n",
    "    return (b - a) * (f(a) + 4 * f((b + a) / 2) + f(b)) / 6\n",
    "\n",
    "\n",
    "def C(f, a, b):\n",
    "    \"\"\"\n",
    "    柯特斯公式\n",
    "    :param f: 被积函数\n",
    "    :param a: 积分下限\n",
    "    :param b: 积分上限\n",
    "    :return: 积分结果\n",
    "    \"\"\"\n",
    "    h = (b - a) / 4\n",
    "    return (b - a) * (7 * f(a) + 32 * f(a + h) + 12 * f(a + 2 * h) + 32 * f(a + 3 * h) + 7 * f(b)) / 90\n",
    "\n",
    "\n",
    "print(T(lambda x: 1 / (x ** 3 + 1), 0, 1))\n",
    "print(C(lambda x: 1 / (x ** 3 + 1), 0, 1))\n",
    "print(S(lambda x: 1 / (x ** 3 + 1), 0, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 复化机械求积公式"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8354045408249832\n",
      "0.8356488543249287\n",
      "0.8356495142270214\n"
     ]
    }
   ],
   "source": [
    "def T_step(f, a, b, n):\n",
    "    \"\"\"\n",
    "    复化梯形公式\n",
    "    :param f: 被积函数\n",
    "    :param a: 积分下限\n",
    "    :param b: 积分上限\n",
    "    :param n: 步数\n",
    "    :return: 积分结果\n",
    "    \"\"\"\n",
    "    h = (b - a) / n\n",
    "    ret = 0\n",
    "    for i in range(1, 1 + n):\n",
    "        ret += 2 * f(a + i * h)\n",
    "    ret -= 2 * f(a + n * h)\n",
    "\n",
    "    return 0.5 * h * (f(b) + f(a) + ret)\n",
    "\n",
    "\n",
    "def S_step(f, a, b, n):\n",
    "    \"\"\"\n",
    "    复化辛普森公式\n",
    "    :param f: 被积函数\n",
    "    :param a: 积分下限\n",
    "    :param b: 积分上限\n",
    "    :param n: 步数 n % 2 == 0\n",
    "    :return: 积分结果\n",
    "    \"\"\"\n",
    "    m = n // 2\n",
    "    h = (b - a) / n\n",
    "    ret = 0\n",
    "    for i in range(1, 1 + m):\n",
    "        ret += 4 * f(a + (2 * i - 1) * h)\n",
    "        ret += 2 * f(a + 2 * i * h)\n",
    "    ret -= 2 * f(a + 2 * m * h)\n",
    "\n",
    "    return h * (f(b) + f(a) + ret) / 3\n",
    "\n",
    "\n",
    "def C_step(f, a, b, n):\n",
    "    \"\"\"\n",
    "    复化柯特斯公式\n",
    "    :param f: 被积函数\n",
    "    :param a: 积分下限\n",
    "    :param b: 积分上限\n",
    "    :param n: 步数 n % 4 == 0\n",
    "    :return: 积分结果\n",
    "    \"\"\"\n",
    "    m = n // 4\n",
    "    h = (b - a) / n\n",
    "    ret = 0\n",
    "    for i in range(1, 1 + m):\n",
    "        ret += 32 * f(a + (4 * i - 3) * h)\n",
    "        ret += 12 * f(a + (4 * i - 2) * h)\n",
    "        ret += 32 * f(a + (4 * i - 1) * h)\n",
    "        ret += 14 * f(a + 4 * i * h)\n",
    "    ret -= 14 * f(a + 4 * m * h)\n",
    "\n",
    "    return 2 * h * (7 * f(b) + 7 * f(a) + ret) / 45\n",
    "\n",
    "\n",
    "print(T_step(lambda x: 1 / (x ** 3 + 1), 0, 1, 16))\n",
    "print(C_step(lambda x: 1 / (x ** 3 + 1), 0, 1, 16))\n",
    "print(S_step(lambda x: 1 / (x ** 3 + 1), 0, 1, 16))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 变步长梯形公式"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8356486098459829\n"
     ]
    }
   ],
   "source": [
    "def T_n(f, a, b, tol=1e-6):\n",
    "    \"\"\"\n",
    "    变步长梯形公式\n",
    "    :param f: 被积函数\n",
    "    :param a: 积分下限\n",
    "    :param b: 积分上限\n",
    "    :param tol: 允许的误差限\n",
    "    :return: 积分结果\n",
    "    \"\"\"\n",
    "\n",
    "    def cal_t():\n",
    "        n = 1\n",
    "        t_n = (b - a) * (f(b) + f(a)) / 2\n",
    "        t_2n = 0\n",
    "        while True:\n",
    "            h = (b - a) / (n * 2)\n",
    "            t_2n = t_n / 2\n",
    "            for i in range(1, n * 2, 2):\n",
    "                t_2n += h * f(a + h * i)  # 注意这里的系数 2\n",
    "            yield t_n\n",
    "            t_n = t_2n\n",
    "            n *= 2\n",
    "\n",
    "    t = cal_t()\n",
    "    t1, t2 = next(t), 0\n",
    "    while True:\n",
    "        t2 = next(t)\n",
    "        if np.abs(t2 - t1) <= tol:\n",
    "            return t2\n",
    "        t1 = t2\n",
    "\n",
    "\n",
    "print(T_n(lambda x: 1 / (x ** 3 + 1), 0, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Romberg 公式"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.8356488482957399\n"
     ]
    }
   ],
   "source": [
    "def Romberg_quadrature_algorithm(f, a, b, tol=1e-6):\n",
    "    \"\"\"\n",
    "    Romberg 公式\n",
    "    :param f: 被积函数\n",
    "    :param a: 积分下限\n",
    "    :param b: 积分上限\n",
    "    :param tol: 允许的误差限\n",
    "    :return: 积分结果\n",
    "    \"\"\"\n",
    "    def cal_r():\n",
    "        def cal_c():\n",
    "            def cal_s():\n",
    "                def cal_t():\n",
    "                    n = 1\n",
    "                    t_n = (b - a) * (f(b) + f(a)) / 2\n",
    "                    t_2n = 0\n",
    "                    while True:\n",
    "                        h = (b - a) / (n * 2)\n",
    "                        t_2n = t_n / 2\n",
    "                        for i in range(1, n * 2, 2):\n",
    "                            t_2n += h * f(a + h * i)  # 注意这里的系数 2\n",
    "                        yield t_n\n",
    "                        t_n = t_2n\n",
    "                        n *= 2\n",
    "                t = cal_t()\n",
    "                while True:\n",
    "                    t1, t2 = next(t), next(t)\n",
    "                    yield 4 * t2 / 3 - t1 / 3\n",
    "            s = cal_s()\n",
    "            s1, s2 = next(s), 0\n",
    "            while True:\n",
    "                s2 = next(s)\n",
    "                yield 16 * s2 / 15 - s1 / 15\n",
    "                s1 = s2\n",
    "        c = cal_c()\n",
    "        c1, c2 = next(c), 0\n",
    "        while True:\n",
    "            c2 = next(c)\n",
    "            yield 64 * c2 / 63 - c1 / 63\n",
    "            c1 = c2\n",
    "    r = cal_r()  # 注意一次循环只用 next() 一次\n",
    "    r1, r2 = next(r), 0\n",
    "    while True:\n",
    "        r2 = next(r)\n",
    "        if (np.abs(r1 - r2) <= tol):\n",
    "            return r2\n",
    "        r1 = r2\n",
    "\n",
    "\n",
    "print(Romberg_quadrature_algorithm(lambda x: 1 / (x ** 3 + 1), 0, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 常微分方程的数值解法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 欧拉法\n",
    "+ 显式欧拉法\n",
    "+ 隐式欧拉法\n",
    "+ 二步欧拉法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1, 1.2, 1.48, 1.8559999999999999, 2.3472, 2.9766399999999997, 3.7719679999999998]\n",
      "[1, 1.2, 1.56, 1.984, 2.5936000000000003, 3.3414400000000004, 4.330176000000001]\n"
     ]
    }
   ],
   "source": [
    "def euler1(f, y0, xs):\n",
    "    \"\"\"\n",
    "    显式欧拉法\n",
    "    :param f: f(x_n, y_n)\n",
    "    :param y0: 初值条件\n",
    "    :param xs: 求值节点，等差\n",
    "    :return: 节点的值\n",
    "    \"\"\"\n",
    "    h = xs[1] - xs[0]\n",
    "    ys = [y0]\n",
    "    for i in range(len(xs)):\n",
    "        ys.append(ys[i] + h * f(xs[i], ys[i]))\n",
    "    return ys\n",
    "\n",
    "\n",
    "def euler2(f, y0, y1, xs):\n",
    "    \"\"\"\n",
    "    二步欧拉法\n",
    "    :param f: f(x_n, y_n)\n",
    "    :param y0, y1: 初值条件\n",
    "    :param xs: 求值节点，等差\n",
    "    :return: 节点的值\n",
    "    \"\"\"\n",
    "    h = xs[1] - xs[0]\n",
    "    ys = [y0, y1]\n",
    "    for i in range(1, len(xs)):\n",
    "        ys.append(ys[i - 1] + 2 * h * f(xs[i], ys[i]))\n",
    "    return ys\n",
    "\n",
    "\n",
    "xs = np.arange(0, 1.2, 0.2)\n",
    "print(euler1(lambda x, y: y + x, 1, xs))\n",
    "print(euler2(lambda x, y: y + x, 1, 1.2, xs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 改进欧拉法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1, 1.24, 1.5768, 2.031696, 2.6306691200000003, 3.4054163264000006, 4.394607918208001]\n"
     ]
    }
   ],
   "source": [
    "def improved_euler(f, y0, xs):\n",
    "    \"\"\"\n",
    "    改进欧拉法\n",
    "    :param f: f(x_n, y_n)\n",
    "    :param y0: 初值条件\n",
    "    :param xs: 求值节点，等差\n",
    "    :return: 节点的值\n",
    "    \"\"\"\n",
    "    h = xs[1] - xs[0]\n",
    "    ys = [y0]\n",
    "    for i in range(len(xs)):\n",
    "        k1 = h * f(xs[i], ys[i])\n",
    "        k2 = h * f(xs[i] + h, ys[i] + k1)\n",
    "        ys.append(ys[i] + (k1 + k2) / 2)\n",
    "\n",
    "    return ys\n",
    "\n",
    "\n",
    "print(improved_euler(lambda x, y: y + x, 1, xs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 龙格-库塔法\n",
    "+ 二阶（中点法、变形欧拉法）\n",
    "+ 三阶（Kutta 公式）\n",
    "+ 四阶（经典龙格-库塔法）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1, 1.24, 1.5768, 2.031696, 2.6306691200000003, 3.4054163264000006, 4.394607918208001]\n",
      "[1, 1.2426666666666666, 1.5833102222222222, 2.043616218074074, 2.6500699410078026, 3.4350187546175297, 4.437969572306209]\n",
      "[1, 1.2428, 1.5836359199999999, 2.0442129126879998, 2.651041651557123, 3.43650227321187, 4.440143876500978]\n"
     ]
    }
   ],
   "source": [
    "def runge_kutta2(f, y0, xs):\n",
    "    \"\"\"\n",
    "    中点法 / 变形欧拉法\n",
    "    :param f: f(x_n, y_n)\n",
    "    :param y0: 初值条件\n",
    "    :param xs: 求值节点，等差\n",
    "    :return: 节点的值\n",
    "    \"\"\"\n",
    "    h = xs[1] - xs[0]\n",
    "    ys = [y0]\n",
    "    for i in range(len(xs)):\n",
    "        k1 = h * f(xs[i], ys[i])\n",
    "        k2 = h * f(h / 2 + xs[i], k1 / 2 + ys[i])\n",
    "        ys.append(ys[i] + k2)\n",
    "\n",
    "    return ys\n",
    "\n",
    "\n",
    "def runge_kutta3(f, y0, xs):\n",
    "    \"\"\"\n",
    "    Kutta 公式\n",
    "    :param f: f(x_n, y_n)\n",
    "    :param y0: 初值条件\n",
    "    :param xs: 求值节点，等差\n",
    "    :return: 节点的值\n",
    "    \"\"\"\n",
    "    h = xs[1] - xs[0]\n",
    "    ys = [y0]\n",
    "    for i in range(len(xs)):\n",
    "        k1 = h * f(xs[i], ys[i])\n",
    "        k2 = h * f(h / 2 + xs[i], k1 / 2 + ys[i])\n",
    "        k3 = h * f(h + xs[i], 2 * k2 - k1 + ys[i])\n",
    "        ys.append(ys[i] + (k1 + 4 * k2 + k3) / 6)\n",
    "\n",
    "    return ys\n",
    "\n",
    "\n",
    "def runge_kutta4(f, y0, xs):\n",
    "    \"\"\"\n",
    "    经典龙格-库塔公式\n",
    "    :param f: f(x_n, y_n)\n",
    "    :param y0: 初值条件\n",
    "    :param xs: 求值节点，等差\n",
    "    :return: 节点的值\n",
    "    \"\"\"\n",
    "    h = xs[1] - xs[0]\n",
    "    ys = [y0]\n",
    "    for i in range(len(xs)):\n",
    "        k1 = h * f(xs[i], ys[i])\n",
    "        k2 = h * f(h / 2 + xs[i], k1 / 2 + ys[i])\n",
    "        k3 = h * f(h / 2 + xs[i], k2 / 2 + ys[i])\n",
    "        k4 = h * f(h + xs[i], k3 + ys[i])\n",
    "        ys.append(ys[i] + (k1 + 2 * k2 + 2 * k3 + k4) / 6)\n",
    "\n",
    "    return ys\n",
    "\n",
    "\n",
    "print(runge_kutta2(lambda x, y: y + x, 1, xs))\n",
    "print(runge_kutta3(lambda x, y: y + x, 1, xs))\n",
    "print(runge_kutta4(lambda x, y: y + x, 1, xs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 方程求根"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 二分法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.3247179985046387\n"
     ]
    }
   ],
   "source": [
    "def bisearch(f, a, b, tol=1e-6):\n",
    "    \"\"\"\n",
    "    二分法求根\n",
    "    :param f: f(x)=0\n",
    "    :param a, b: 近似根区间\n",
    "    :tol: 允许的误差限\n",
    "    :return: 方程的近似根\n",
    "    \"\"\"\n",
    "    a, b = min(a, b), max(a, b)\n",
    "    while tol < b - a:\n",
    "        m = (b + a) / 2\n",
    "        if f(a) * f(m) < 0:\n",
    "            b = m\n",
    "        else:\n",
    "            a = m\n",
    "\n",
    "    return (b + a) / 2\n",
    "\n",
    "\n",
    "print(bisearch(lambda x: x ** 3 - x - 1, 1, 1.5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 不动点迭代法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.3247177357599116\n"
     ]
    }
   ],
   "source": [
    "def iter_method(g, x0, tol=1e-6):\n",
    "    \"\"\"\n",
    "    不动点迭代法方程求根\n",
    "    :param g: x=g(x)\n",
    "    :param x0: 根的初始值\n",
    "    :tol: 允许的误差限\n",
    "    :return: 方程的近似根\n",
    "    \"\"\"\n",
    "    g1, g2 = x0, 0\n",
    "    while True:\n",
    "        g2 = g(g1)\n",
    "        if tol >= np.abs(g2 - g1):\n",
    "            return g2\n",
    "        g1 = g2\n",
    "\n",
    "\n",
    "print(iter_method(lambda x: (1 + x) ** (1 / 3), 1.3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 牛顿法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.324717957244781\n"
     ]
    }
   ],
   "source": [
    "def newton_method(f, f_diff, x0, tol=1e-6):\n",
    "    \"\"\"\n",
    "    牛顿法方程求根\n",
    "    :param g: x=g(x)\n",
    "    :param x0: 根的初始值\n",
    "    :tol: 允许的误差限\n",
    "    :return: 方程的近似根\n",
    "    \"\"\"\n",
    "\n",
    "    x1, x2 = x0, 0\n",
    "    while True:\n",
    "        x2 = x1 - f(x1) / f_diff(x1)\n",
    "        if tol >= np.abs(x1 - x2):\n",
    "            return x2\n",
    "        x1 = x2\n",
    "\n",
    "\n",
    "print(newton_method(lambda x: x ** 3 - x - 1, lambda x: 3 * x ** 2 - 1, 1.9))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
